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I. INTRODUCTION 

Spherical, three dimensional dusty plasma crystals 
(Yukawa balls) can be created in laboratory experiments 
by trapping dust particles inside a glass box above the 
lower electrode of an rf discharge [TJ [2] . They can be re- 
garded as the analog to confined Coulomb clusters in ion 
traps [3j H] with the difference being the screened parti- 
cle interaction [5 j. Recent experimental and theoretical 
work has been concerned with their ground state configu- 
rations [6]{9] , the occurrence of metastable states [T0HT2] , 
time dependent shell formation [13] and the investigation 
of their normal modes [14-16 , for an overview see [T7] . 

Collective modes contain information about the par- 
ticle interactions and the external confinement condi- 
tions. They are of similar interest for many-particle sys- 
tems as are the atomic spectra in atomic and molecular 
physics. By comparing with theoretical predictions ex- 
perimental normal mode measurements can be used to 
infer, for example, the confinement parameters and the 
particle charge from experimental data without changing 
the plasma conditions p~4] [18] . Depending on the quality 
of the theory this makes them a valuable non-invasive 
tool for experimentalists. 

In this paper we study the collective excitations of a 
spherically confined Yukawa plasma applying three in- 
dependent methods: molecular dynamics (MD) simu- 
lations, normal mode analysis and a recent continuum 
approach [19 based on cold fluid theory. In Ref. [19] 
it was predicted that a spherically trapped plasma with 
Yukawa interaction possesses a much richer collective ex- 
citation spectrum than its Coulomb interacting counter- 
part. In particular, it was found that there exists an 
additional class of modes characterized by a mode in- 
dex n = 0, 1, 2, . . . corresponding to the number of radial 
nodes in the perturbed electrostatic potential inside the 
plasma. 

The main goal of this paper is to verify these theo- 
retical predictions. The success of the fluid theory for 
Coulomb systems has been demonstrated both theoret- 



ically [20] and experimentally [4] [21]. However, recent 
work on the ground state density profile of a Yukawa 
plasma [22] [23], on which the fluid theory is based, indi- 
cates that correlation effects tend to become more impor- 
tant in systems where the particle interaction is screened 
and the interaction becomes short-ranged. Thus, a sys- 
tematic investigation and detailed comparison with exact 
simulations is required in order to explore the applicabil- 
ity limits of the fluid theory also for Yukawa plasmas. 

By comparison with first principle MD simulations 
we prove that the additional modes are physical real- 
ity. Based on the very good agreement between MD and 
experiment [SJ [12], we expect that they should be ob- 
servable in Yukawa ball experiments as well. Performing 
detailed MD simulations covering a broad range of cou- 
pling and screening parameters we are able to determine 
the applicability range of the fluid theory: It predicts 
the frequencies of the low order modes with n = with 
good accuracy and those of the n = 1 modes for weak 
to moderate screening. For n > 2 or strong screening 
the deviations quickly become larger. The main source 
of the observed deviations are correlation effects which 
are missing in the cold fluid theory. 

Among the eigenmodes especially the breathing mode 
has attracted considerable attention in recent years [15] 
[24] [25] . We therefore also compare the fluid breathing 
mode (the lowest monopole mode) with its counterpart 
among the eigenmodes of the discrete iV-particle system 
in the crystalline state and the results of MD simulations. 

This paper is organized as follows. In Sec.[ll]we briefly 
recall the main results of the fluid theory [19] which are 
relevant for the present work. The simulation method is 
explained in some detail in Sec. |III| an d the results are 
compared with the theory in Sec. |IV| Section [V] deals 
with the exact eigenmodes of the crystallized iV-particle 
system and examines the relations of the breathing mode 
found in cold fluid theory with the crystal eigenmodes 
and the MD simulations. We conclude with a summary 
of our results in Sec. IVI1 
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II. REVIEW OF FLUID MODES 

In this section we briefly review recent results for the 
ground state and normal modes of a Yukawa plasma in 
an external confinement V(r) = muj\r 2 /2 to make this 
paper as self-contained as possible. 

For Coulomb interaction the ground state density pro- 
file in the cold fluid approximation is a sphere of con- 
stant density while for Yukawa interaction the density 
exhibits a parabolic decay towards the plasma boundary 
R(N, k) [22] . Here n denotes the inverse screening length. 
The radius must be determined from the following equa- 
tion for £ = kR, 



£ 6 + 6£ 5 + 15 [£ 4 + £ 3 -fc3(£ + l)] =0) 



(1) 



where kc = kRq and Rq = aN 1 / 3 denotes the plasma 
radius for a Coulomb system. The relevant length scale 
is a = (q 2 /muj 2 ) 1 ^ , corresponding to the Wigner-Seitz 
radius in the Coulomb limit. Accurate analytical approx- 
imations for £ are available [19]. 

The normal modes have been derived in Ref. [19] by 
linearizing the cold fluid equations. As a result the first 
order potential perturbation is obtained as 



^ n/out (r,u,) = / in /° ut (r,^r(^), 



(2) 



where the indices in/out denote the regions inside (r < 
R) and outside (r > R) the plasma. 

The radial function inside the plasma is given by 



a£ - 8z a£ + « 



2 2 



, (3) 



where the parameters of the hypergeometric function 2^1 
are 



<*i = e + ^, ^ = y^ + i) + ^ + 2^ 2 , 

and x 2 = 2(^(0) — ft 2 ). We use the notation ft = uj/ujq 
for the scaled eigenfrequency. Further, x — kt and 



uj 2 p {r) 



C 2 3 + g x 2 
2 l + £ 2 



denotes the local plasma frequency, which is determined 
by the ground state density profile. 

Outside the plasma the radial function takes the form 
of a modified spherical Bessel function, 



/ out (r) ~ k e ( K r). 



(4) 



The eigenfrequencies are determined by matching the 
two solutions at the plasma boundary. They satisfy 
ft 2 < i.e. the plasma frequency at the boundary 

is the maximum frequency for regular normal mode solu- 
tions with a discrete spectrum. For ft 2 (0) > ft 2 > ft 2 (^) 
the spectrum should become continuous, for details cf. 
Ref. [26]. 

Let us now summarize the properties of the normal 
mode spectrum which follows from cold fluid theory [19] : 



1. Normal modes are characterized by three mode 
numbers: the radial mode number n = 0, 1, 2, . . . 
and the two angular mode numbers £ and m with 
\m\ < £. Their frequencies solely depend on n and 
the dimensionless plasma parameter £, see Eq. ([!]). 

2. The n = modes are directly related to the surface 
modes of the associated Coulomb system [20l [271 
[28] . where ft 2 = 31/ (2£ + 1), while the remaining 
modes with n > 1 originate from its degenerate 
bulk modes with ft 2 = 3. Eigenmodes with n = 
and arbitrary £ and n = 1 with £ = 0, 1, 2 exist for 
all £ whereas all other modes exist only for £ > £ crit 



where 



£-crit 



1 

2 

(2n 



Cn£ + V CniiCni + 4) 

l){n + l) -4. 



nt ' 



(5) 



At these points the plasma frequency assumes the 
values 



^(^f) = (2n+l)(n-l) + (2n-l)£ 



(6) 



3. In the limit £ —> oo, corresponding to strong screen- 
ing or the macroscopic limit with N — >• oo, the 
mode frequencies are known analytically and read 



O 2 



2n 2 + (2^ + 3)n + l 



(7) 



For finite £ they must be determined numerically 
from Eq. (29) in Ref. [19]. 



III. MOLECULAR DYNAMICS SIMULATIONS 

A. Simulation Method 

A molecular dynamics (MD) simulation allows for an 
investigation of the normal modes without approxima- 
tions with respect to the particle interaction or their 
temperature and will serve as a benchmark for the fluid 
theory. 

In order to simulate the same conditions as in the fluid 
approach, we consider N particles interacting through 
a Yukawa pair potential <j>{r) = q 2 e~ Kr /r in an exter- 



nal confinement V(r) — 
Hamiltonian is given by 



qT /2. The corresponding 



N 2 N i N 

^Ei^ + E^D + ^E^D- (8) 



U (ri,...,rjv) 



In dusty plasma experiments collisions with neutral par- 
ticles can be an important factor for the dynamics. We, 
therefore, include an additional damping term and a fluc- 
tuating force (Langevin dynamics) in the equation of mo- 
tion, which is given by 



mri = — ViE/(ri, . . . ,r N ) - mvvi + f^(t). 



(9) 



3 



In thermodynamic equilibrium, the friction coeffi- 
cient v and the Gaussian noise £j(t) are related by 

[if (t)ff(t')) = ZmvkBTSijSapSit - t f ), where G 
{x,y,z} and z,j G {l,...,iV}. The system is fully 
characterized by the particle number A/", the (normal- 
ized) screening parameter kcl, and the coupling parame- 
ter T = q 2 /(ak B T). 

Note that the so-defined V is not the same as in a 
macroscopic, homogeneous system because of the inho- 
mogeneous density in the trap. Only in the Coulomb 
system the mean density, and thus the Wigner Seitz ra- 
dius, is constant. For Yukawa systems V can be regarded 
as a dimensionless inverse temperature rather than a cou- 
pling parameter. Here the effective coupling varies, de- 
pending on the distance from the center of the trap and 
the screening parameter. Effective coupling parameters 
which take screening into account have been analyzed for 
two- [29 , 30] and three-dimensional [31] Yukawa systems. 

We perform finite temperature simulations where the 
normal modes are thermally excited and no external ex- 
citation or driving mechanism is required. When the 
Langevin method is used (y ^ 0), equilibration of the 
system is ensured by friction and the stochastic forces. 

For v = we rescale the individual particle velocities 
in the equilibration phase to enforce the desired coupling 
parameter via the relation (E^ n ) = 3iV^T/2. Addi- 
tionally, we remove any residual rotation and center of 
mass motion of the whole plasma. After this initial equi- 
libration period the system is evolved microcanonically. 



B. Mode detection 

In order to make contact with the fluid modes we need 
to find a way of measuring the fluid oscillations in the MD 
simulation. Since the fluid modes ([2| are perturbations of 
the induced potential, we have to derive an analogous ex- 
pression which is appropriate for the MD simulation. The 
same method was used in [20 for a spheroidal Coulomb 
system. 

For a general charge distribution qn(r, t) the solution 
of the (screened) Poisson equation [22] 

(A - K 2 )^ tot = -Airqn, (10) 

with the boundary condition tot (r,£) — >> at infinity is 

-k\y— y' I 



"tot 



(r,t)=q 



f e~ K 



-dr'. 



(11) 



We can now expand the Yukawa potential, which is the 
Green's function of Eq. (10) for this boundary condition, 



in a series of spherical harmonics according to [32] 

6 , [ ,.' = 4tt« ^(w<)A*(«r>) y/ m (0', <t>')Y lm {e, (/>). 

' r r ' im 

The expansion involves the modified spherical Bessel 
functions it(x) and ki(x) with the familiar arguments 
r< = min(|r|, |r'|) and r> = max(|r|, | r / 1 ) . 



With this expression at hand we can define multipole 
moments of the density. This is easily achieved by using 
the Green's function expansion in Eq. ( 11 ) with the result 



Qim(t) 



47T 



21- 



n(r',t) i^r')Y[ m (6' ^')dv' . 



Here we have introduced 



(2^+1)!! 



(2^-1)!! 



(12) 



where the prefactors are chosen such that the multipole 
moments reduce to their usual form in the limit k —> 0. 

These definitions allow us to write the general solution 
for the potential, Eq. ( [TT] ), outside the region where the 
density is non-zero (we used r< as integration variable) 
as 



t(r,t) = £ 



im 



47T 
2ZTl 



qem(t) ke(K, r) Y irn (0 : cj)). (13) 



The full time-dependence of 0t o t(r, i) is now contained 
in the multipole moments ( 12 ). They are easily accessible 



in the MD simulation because the integrals in Eq. (12) 



reduce to a sum over all particles. They contain the re- 
quired information about the oscillations of the potential 
with the mode numbers £ and m. 

The mode frequencies are extracted from qe m (t) by 
computing the associated power spectral density 



Qlm{u) = lim 

T— >-oo 



\F(qe m (t)/N,T,w)f 



(14) 



where F (q£ m (t) / N , T, uj) denotes the Fourier transform 
of q£m(t)/N over a finite time interval T [38] . Maxima 
of Qim(u) indicate frequencies at which the system sup- 
ports collective excitations. 



IV. MD SIMULATION RESULTS 

When comparing the MD simulation and the fluid the- 
ory we immediately recognize several differences. The 
MD simulation is performed at finite temperature while 
the fluid modes correspond to a system at T = 0, i.e. in- 
finite coupling strength. Furthermore, the fluid approach 
neglects correlation effects which are fully included in the 
simulation. These limitations must be kept in mind when 
comparisons between the two are made. 

The spectra obtained from the MD simulations confirm 
the results of cold fluid theory, e.g. the existence of addi- 
tional modes, but the quantitative agreement depends on 
several parameters and the particular mode. In the fol- 
lowing we will, therefore, separately discuss the influence 
of the particle number, screening, coupling and friction 
on the mode spectra and compare the simulations with 
the fluid approach. 
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FIG. 1: Screening dependence of the m = multipole spectra 
Qm(u) (arb. units) for various t with A/" = 1000 (left column) 
and N = 4000 particles (right column) at F = 150 without 
friction. The solid lines show the results of cold fluid theory. 



1. Screening dependence 

We begin our discussion with the influence of the 
screening parameter on the normal modes. Figure [I] 
shows mode spectra for N = 1000 and 4000 particles at 
strong coupling conditions and their dependence on kcl. 
Note that for the same T the effective coupling of the sys- 
tem can be quite different, depending on kcl, e.g. [501155]. 
On the one hand, the interaction strength decreases as 
the screening parameter is increased, leading to a weaker 
effective coupling. On the other hand, as the interaction 
range 1/n is reduced, the density increases at the same 
time. The cluster is being compressed which, in turn, in- 
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FIG. 2: a) Monopole and b) dipole spectra for various screen- 
ing parameters with N — 1000 particles at Y = 150. The 
thick arrow in b) indicates the low frequency peak in the 
dipole spectrum. 



creases the effective coupling. Thus, these two effects act 
in opposite directions and may at least partially cancel. 

In the Coulomb limit, na <C 1, we observe max- 
ima of Q^m(cj) close to the surface mode frequencies 
Qi = <\/3£/ (2£ + 1). The center of mass motion was re- 
moved in the equilibration phase so there is no maximum 
at fii = 1 for £ = 1, see Fig.[l]3. When friction is included 
we observe a strong peak at the expected position. Ad- 
ditionally, there are maxima near the plasma frequency 
ft = y/S, except for £ = 0, where the multipole moment 
is constant (for n = the multipole moment Qoo is smi_ 
ply the total charge of the cluster). The Coulomb modes 
have been investigated in detail in [20] . For Yukawa inter- 
action the general trend of our simulation results is that 
all mode frequencies increase with na and new modes ap- 
pear at strong screening which is the expected behavior 
from cold fluid theory. Let us now compare the results 
in more detail. 

Monopole modes. For £ = [see Figs. [T^i,e] there is 
no mode with n = 0. Comparing the results of cold fluid 
theory with the simulations we find excellent agreement 
for the lowest mode with n = 1. There are only small de- 
viations which increase with kcl. The second mode with 
n — 2 exists only for £ > £30* ~ 2.73 in cold fluid theory, 
corresponding to na ~ 0.35 for N — 1000 and na w 0.22 
for N = 4000. This result is confirmed by the MD sim- 
ulation. However, there is a small upshift of the MD 
frequencies compared to the fluid theory which is greater 
for N = 1000 than for N = 4000. For both particle num- 
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bers we observe a third mode in the simulations at even 
higher frequencies which should correspond to the n = 3 
mode of cold fluid theory. However, for N = 1000 this 
mode does not yet exist in the theory at the simulated 
screening parameters and only appears for na > 2.88. 

Since the fluid theory is a continuum approach it is ex- 
pected to yield accurate results only for large N and the 
observed discrepancies are easily understood. In particu- 
lar, cold fluid theory should yield reliable results only for 
modes that oscillate on large length scales, i.e. for long 
wavelength modes. In the present case with £ = m = 
the length scale of the oscillations is mainly determined 
by the number of radial nodes, i.e. the mode number n. 
While for n = 1 the agreement between simulation and 
theory is very good, the deviations rapidly increase as 
we go to n > 2, due to the short wavelength character 
of these modes. This is qualitatively similar to the situ- 
ation in a macroscopic Yukawa plasma where the Vlasov 
(mean field) model for the longitudinal dispersion rela- 
tion Lj(k) in the strongly coupled phase is only adequate 
for ka <C 1, where k is the wave number [34]. Here 
the inclusion of correlation effects via the quasi-localized 
charge approximation was shown to be essential for a reli- 
able description of the dispersion relation for larger wave 
numbers. 

Dipole modes. For £ = 1 the center of mass mode 
is an exact eigenmode with ft = 1, independent of na. 
Since this result is confirmed by the theory this special 
mode is only of minor interest here. More interesting are 
the higher order modes with n = 1 and 2 [see Figs. [lj>,f]. 
Similar to the results for £ = we observe good agree- 
ment for n = 1. Again, deviations between the two ap- 
proaches increase with screening. The general behavior 
for n = 2 is the same. The mode exists in the MD sim- 
ulation only at strong screening and the starting point 
is well captured by the theory. Deviations between the 
mode frequencies are larger than for £ = and, again, 
smaller for N = 4000 than for TV = 1000. A striking dif- 
ference compared with £ = is the existence of another 
low frequency mode which has no analog in cold fluid 
theory. Its frequency is almost independent of na. For 
N = 4000 the frequency is lower than for 1000 particles. 

Quadrupole and octupole modes. In the case of 
£ = 2 and £ — 3 we find the first non-trivial modes with 
n = 0. Here the agreement between simulation and the- 
ory is excellent. As in the previous cases the (small) 
deviations increase with na. This also holds for n = 1. 
For n = 2 the deviations between the theoretical mode 
frequencies and the simulation become rather large. The 
starting point for the existence of the modes is still well 
reproduced by fluid theory, at least for N = 4000. Again, 
there are low frequency excitations below the n = 
modes which do not exist in cold fluid theory. 

Inspecting the mode spectra more closely we observe 
another important property of the modes that exist only 
beyond a certain screening parameter in cold fluid theory. 
Left to their critical points these modes reappear in the 
simulations at lower na for a finite screening interval, but 



with a relatively low intensity. This is shown in Fig. [2] in 
more detail. In the monopole spectrum [Fig. [2^t] there is 
a third maximum at Q ~ 4.34 for na — 2 which is absent 
for na — 1.4 but reappears at Q « 2.99 for na — 1. 
Fig. [2}d shows the same behavior for the dipole modes. 
Here the second peak at Q « 3.18 for na = 1 is absent 
in the na = 0.7 spectrum. For na = 0.5 it reappears at 
Q « 2.6 as a second maximum rather close to the main 
peak. Since cold fluid theory does not allow for additional 
normal mode solutions below fi p (£) these excitations may 
be a due to the continuous part of the spectrum. 

Another possible explanation could be the density pro- 
file in the MD simulation which differs quite substantially 
from the cold fluid result. While at the large coupling 
parameters we used in the simulation the exact density 
has a shell structure, cold fluid theory only describes 
the mean density correctly. Improvements in this regard 
are possible by including correlations via the hypernetted 
chain approximation [35j[36|, see Sec. 



VI 



2. Influence of coupling strength 

Next, we consider temperature effects on the mode 
spectra. Fig. [3] shows the same spectra as Fig. [I] but for 
T = 30. Even though the coupling is five times weaker 
in this case the particles are still strongly coupled. The 
main difference is a structural change of the plasma. For 
T = 150 the density profile shows a clear shell structure 
while for T = 30 the density profile is well described by 
the cold fluid mean field result, except for the boundary 
region where it smoothly goes to zero [22] . Additionally, 
the density has small maxima near r = R which are first 
indications for the onset of the shell structure [T3l [28] . 
This is a correlation effect not included in cold fluid the- 
ory [35| [36] . 

Comparing the mode frequencies with those for T = 
150 we find that, on the one hand, they are only weakly 
affected by the change of T. We still observe good agree- 
ment with cold fluid theory as in the previous case. At 
lower coupling the peak width increases which is indica- 
tive of higher damping. On the other hand, there are 
two important qualitative changes of the spectra. First, 
the low frequency modes for £ — 1,2,3 disappear com- 
pletely. Second, the modes that exist only for ^ n £ > 
do not reappear at lower screening. For the monopole 
and dipole modes this can be seen more clearly by com- 
paring Figs. [2] and [4] which show the same spectra at 
different coupling strengths. Clearly, the low frequency 
peak in the dipole spectrum is absent at T = 30 (cf. thick 
arrow in Figs. J2|d, [4)3). Further, the higher order modes 
with n > 1 are strongly damped. 

Figs. [5] and [6] show the coupling dependence in more 
detail. Fig. [5] confirms that the frequencies are almost 
independent of the coupling parameter, at least in the 
strong coupling regime. While for T > 100 the peaks 
are well defined, for T < 100 the peaks become broader 
and less pronounced. The vanishing of several modes is 
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FIG. 3: Screening dependence of the m = multipole spectra 
Qm(u) (arb. units) for various £ with N = 1000 (left column) 
and N = 4000 particles (right column) at F = 30 without 
friction. The solid lines show the results of cold fluid theory. 



clearly visible in Fig.[6j At very strong coupling, F = 400, 
the spectrum shows five distinct peaks. Decreasing T the 
peaks become broader and the higher order modes van- 
ish. The rather broad maximum at low frequencies is 
still visible at T = 121 but is absent in the T = 25 spec- 
trum. The peak positions are only weakly affected by the 
variation of T. While the frequency of the n = mode 
slightly decreases as T decreases, the mode frequency for 
n = 1 increases for £ — 3 and decreases for £ = 2. We did 
not investigate this dependence for all modes in further 
detail because the changes are small and our main inter- 
est is in the strongly coupled region with T > 100. More 
details on the T dependence of the breathing mode can 




FIG. 4: Same as Fig.|2]for F = 30. The thick arrow indicates 
the position of the low frequency mode for F = 150. 



be found in Sec. El 



3. Influence of damping 



So far our results have been obtained for zero damp- 
ing. In order to make predictions for experiments a finite 
damping term must be taken into account. A typical ex- 
ample is shown in Fig. [7| For v = the spectrum has 
two main peaks at Q = 1.32, 2.48 and two smaller peaks 
at Q « 0.3, 2.8. With the increase of v the peak width 
increases and their height decreases. At high frequencies 
we see an increase of the noise level. 

For v/ujq = 0.001 all four peaks are still present in the 
spectrum. At u/ojo = 0.1 the two smaller peaks cannot 
be identified any more and at v/ujo = 1 the spectrum 
is almost fiat and only the main peak can be observed. 
This shows that a low damping rate is essential in or- 
der to observe the higher modes with n = 2, 3, . . . in a 
Yukawa system. In dusty plasma experiments the fric- 
tion coefficient can be within a wide range [TDJ [14j [37] , 
depending on the gas pressure. We thus expect that it 
should be possible to observe the modes in experiments. 
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FIG. 5: Influence of the coupling strength on the m = 
multipole spectra Qio (arb. units) for N — 1000 and 10 < 
T < 400. The screening parameters are kcl = 0.6 (a,b) and 
na — 2 (c,d). 




FIG. 6: Quadrupole and octupole spectra for N — 1000 and 
kol — 2 for three values of the coupling parameter Y as indi- 
cated in the figure. 



V. EXACT CRYSTAL NORMAL MODES AND 
COMPARISON WITH FLUID THEORY AND MD 
SIMULATION 

The collective modes of fluid theory arise from a con- 
tinuum model of the plasma which completely neglects 
the underlying discrete particle structure. On the other 
hand, N strongly coupled correlated particles in a trap 
possess 37V normal modes (related to phonons) which can 
be directly computed applying harmonic lattice theory. 
It is, therefore, very interesting to compare the spectra 
following from fluid theory to the normal mode spectrum. 

A. Frequency sum rule 

Before comparing the two approaches (fluid theory 
and exact normal modes) we discuss a sum rule for 
the squared eigenfrequencies which is known to be par- 
ticularly simple in Coulomb systems. The exact nor- 
mal modes of the discrete TV-particle system are calcu- 
lated from the Hessian matrix U s of a stationary state s 
(ground or metastable state). The total potential energy 
is 

N N 

^ = E^o) + 2E^I4-,ol)> (15) 

where the {r| } are the particle positions. The eigen- 
values of U s are the (squared) eigenfrequencies of the 
system, e.g. [T5] . 




uj/ujq 



FIG. 7: Influence of a finite damping coefficient on the 
quadrupole spectrum for the parameters indicated in the fig- 
ure. Results are shown for is/ujo = 0, 0.01, 0.1, 1. 



A frequency sum rule for particles with Coulomb in- 
teraction in a harmonic trap V(r) = muj^r 2 /2 is given 
by Ei5iK s ) 2 = 3iVcjg [20] . Particularly, the sum is in- 
dependent of the stationary state used to evaluate the 
Hessian. We will show below that this is not a general 
rule and does not hold for the Yukawa system considered 
here. 

Since, in general, Y^i=i m { UJ t ) 2 = Tt(U s ) we need to 
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evaluate the diagonal elements of the Hessian, 



N 



Tr(W) = ^A i 



N N 



k=l 



(16) 



r| . For the harmonic confinement considered 



3mco>Q. The interaction term 



at r, 

here we have AV(r) 
is simply the radial part of the Laplacian applied to 
<j)(r). As the Yukawa potential satisfies (A — K 2 )e~ Kr /r 
— 47r5(r), the result is 



3N 



37V + 2 (nof 



q 2 /a 



(17) 




100 200 
mode number j 



300 



400 800 
mode number j 



1200 



FIG. 8: Projection of the n — 1 monopole mode {£ = m = 0) 
of fluid theory on the 3 N crystal eigenmodes. Results are 
shown for na — 0.4, 1.2, 2 (from left to right peak). The modes 
are ordered by their frequency (from high to low). 



where U? nt denotes the total interaction energy of the 



particles in the state s [second term in Eq. (15)]. The 
delta function does not contribute since 



^0. 



For k = Eq. (16) reduces to the known Coulomb 
limit. Due to the additional interaction term the fre- 
quency sum for Yukawa interaction explicitly depends on 
the particle configuration, e.g. the shell occupation num- 
bers or the arrangement of the particles on a given shell. 
It is not a universal quantity like in the Coulomb case. 
Eq. (17) connects the eigenfrequencies of the state s with 



its interaction energy, i.e. it establishes a connection be- 
tween dynamic and static properties. We computed both 
sides of Eq. (17) independently and found that the sum 
rule is indeed fulfilled. 



B. Comparison of MD simulation with crystal 
eigenmodes 

In the crystalline phase the particle motion is mainly 
determined by the crystal eigenmodes and, consequently, 
the same applies to the oscillations of the multipole mo- 
ments. In the following we extract the dominant contri- 
bution to the multipole moments by identifying the eigen- 
mode that is closest to the respective fluid mode. This is 
achieved by computing the particle displacements accord- 
ing to cold fluid theory and projecting all 3N eigenmodes 
on the so-defined reference mode, see also [20 j. This is 
done according to 



w(ri,...,rjy) 
|w(ri,...,rjv)| 



where Vj is the jth (normalized) eigenvector of the sta- 
tionary state of the crystal with particle positions {r| } 
and w is the eigenvector calculated from cold fluid the- 
ory. The particle displacements for the monopole modes 
in this case are Wj = /'(r, u)e r \ r=rS . The full displace- 
ment vector is then taken as w = (wi, . . . , wjv)- The 
mode with the greatest projection is chosen for compari- 
son with the simulation data. The stationary states used 
for the calculation of the normal modes were found by 



slowly cooling a random initial state in a molecular dy- 
namics simulation until the equilibrium positions were 
reached, e.g. [6] [12]. This procedure was repeated 200 
times for each parameter set (TV, na). The state with the 
lowest energy was then chosen to evaluate the Hessian. 

The projections on the monopole mode with n = 1 
are shown in Fig. [8] One can clearly see that only very 
few modes have a non-negligible projection on this fluid 
mode. There is a single mode with a high projection - 
this corresponds to a 'breathing mode'. While for low 
na the breathing mode is among the normal modes with 
the highest frequencies, the peak shifts to intermediate 
frequencies as screening increases. This is very similar to 
the behavior in cold fluid theory where additional modes 
appear at larger screening which have frequencies higher 
than that of the breathing mode, cf. Fig. [I] In contrast, 
for Coulomb interaction the breathing mode is always the 
one with the highest frequency [20] . 

The very narrow peak of the breathing mode in the 
MD simulations allows for a detailed investigation of its 
frequency. In Fig. [9] we compare the frequency obtained 
in the MD simulation with the frequency of the projected 
crystal eigenmode [peak in Fig. [§]. While for low na the 
frequency in the simulation exceeds the eigenfrequency 
of the crystal, the trend is reversed for larger na. As the 
coupling parameter is increased the deviations between 
the two frequencies decrease and vanish almost identi- 
cally in the very strong coupling limit, V — 150. 

The increase of the deviations with temperature (re- 
duction of T) is readily understood by the limitations of 
harmonic lattice theory: when T increases, the amplitude 
of particle fluctuations around their stationary positions 
grows and cannot be described by the second order con- 
tribution in the expansion of the potential energy. An- 
harmonic corrections and mode coupling effects (which 
are contained in the MD simulation) lead to deviations 
from the normal modes of the harmonic approximation. 
The interesting observation following from Fig. [9] is that 
the sign of these deviations changes with na. We note 
that we also compared the breathing frequencies of dif- 
ferent stationary states for a few cases and found that the 
typical difference is of the order AQ ~ 10 -4 and therefore 
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FIG. 9: Deviation of the n = 1 monopole mode from the 
frequency of the exact eigenmode with the greatest projection 
on the fluid mode for a) N = 100 and b) N = 400 at T = 
5, 10, 30, 150. The solid lines show a linear fit. 





0.00 * 






1 — 1 








1 


-0.01 
















-0.02 






3° 




N-- 






-0.03 


100 


f 




400 -^»— 




-0.04 


- 1000 — 



a) 



Ka = 0.6 


b) ■ 


Ka = 


= 1.5 ; 








na = 2 : 







1 

h'd 



1.5 



200 400 600 800 1000 

TV 



FIG. 10: Frequency deviation of the exact eigenmode with 
the greatest projection on the fluid mode from the results of 
cold fluid theory as a function of a) screening and b) number 
of particles. 



about an order of magnitude smaller than those caused 
by the finite temperature in the MD simulation. 

The observed temperature dependence is peculiar: At 
low screening the MD mode frequency decreases with an 
increase of T while it increases at larger na. This behav- 
ior is in qualitative agreement with the results of [25J who 
found that the breathing frequency increases (decreases) 
with the coupling strength for short- (long-) ranged forces 
of the power law type. Even though their analysis does 
not directly apply to our results with Yukawa interac- 
tion the general trend is the same. In our system the 
interaction range is essentially determined by the param- 
eter £ = kR. In a very small interval around na ~ 0.6 
(N = 100) and Ka « 0.4 (N = 400) the breathing fre- 
quency is almost independent of T. Here £ ^ 2.3 and 
£ « 2.4, respectively, which could be an indication that 
the temperature dependence of the breathing frequency 
in the strongly coupled regime is essentially determined 
by the value of £. 



C. Comparison of fluid theory with crystal 
eigenmodes 

The previous analysis showed that the frequency in the 
MD simulation approaches the frequency of the crystal 



eigenmode in the very strong coupling regime. In order 
to test the validity of the fluid theory, which was derived 
for the infinite coupling case (T — >> oo), we, therefore, 
compare its frequencies with those of the exact crystal 
eigenmodes. This has two advantages over the MD sim- 
ulation: 1) the mode frequency can be determined with 
a high accuracy and 2) the infinite coupling condition is 
already 'built-in'. 

The results are shown in Fig. [lOj In the Coulomb 
limit the breathing mode is among the bulk modes with 
the (exact) frequency Q = y/3 and uniform particle 
displacements [15]. Even though the particle displace- 
ments in the n = 1 monopole mode for £ — >• in 
cold fluid theory do not become uniform [13], the mode 
with the greatest projection is nevertheless the breathing 
mode, and so the frequency deviations in the Coulomb 
limit vanish, at least for the particle numbers considered 
here. For finite screening no universal uniform breath- 
ing mode exists [15]. Here we compared the projection 
according to the theory with uniform displacements, i.e. 
Wj = re r L _« , and found that both methods select the 

r— r i,0 

same modes (with an exception for N = 100, na = 0.1). 
In almost all cases the projection in the latter case is 
slightly less than in the former, i.e. the mode form is 
more accurately described by non-uniform particle dis- 
placements. 

The deviations of the mode frequencies are of the or- 
der of a few percent and about an order of magnitude 
larger than those caused by finite temperature effects, 
cf. Figs. [9] and [lOj They increase linearly with na and 
decrease with TV. It was shown in [24] that the breath- 
ing frequency increases if correlations are included in the 
theory. This explains why the fluid treatment underesti- 
mates the mode frequencies. Our results suggest that in 
the limit N — >• oo a crystal eigenmode could exist that 
directly corresponds to the n = 1 monopole fluid mode. 
Its frequency would be given by tt = Vb [cf. Eq. Q for 
(n, i) = (1,0)], in contrast to the Coulomb limit where 

n = Vs. 



VI. CONCLUSION 

In this paper we have analyzed the normal modes 
of a spherically confined Yukawa plasma over a broad 
range of coupling parameters T and screening parame- 
ters kcl. This question was studied using three differ- 
ent approaches: MD simulations, harmonic lattice theory 
(crystal normal modes) and a recently developed fluid 
theory. Of central interest was to test the predictions of 
the latter: the existence of an additional (compared to 
Coulomb systems) class of collective modes enumerated 
by the mode number n. 

Our first conclusion is that 1) these modes are also 
present in the exact results. 2) The lowest order mode 
(n = 0) of fluid theory practically coincides with the 
MD result for all studied cases. 3) With increasing mode 
index n > 1 the frequencies predicted by fluid theory 
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are systematically too low whereas the minimal value kcl 
where the modes exist is too large. 4) The results of fluid 
theory become more accurate when the particle number 
increases, in agreement with the continuum character of 
the approach. 5) Comparison of the fluid monopole mode 
with n = 1 to the breathing mode of crystal harmonic 
lattice theory showed that the main source of error are 
missing correlations: The frequency of fluid theory is sys- 
tematically too small and deviations increase with kcl. 6) 
Stimulated by the predictions of fluid theory we have per- 
formed detailed MD simulations which have uncovered 
the complete collective excitation spectrum. In addition 
to the fluid theory we observed two peculiarities which 
are present at strong coupling: i) a low frequency mode 
below the n = modes of fluid theory and ii) reentrant 
character of fluid modes (which exist beyond a critical 
value of n) at low values of na. The physical origin of 
these effects is still unclear and will be explored else- 
where. 

Since the effects of a finite temperature on the mode 
frequencies is negligible when compared with those 
caused by correlations, an improved fluid theory should 
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